m = 2;
lambda = 0.1;
mu = 0.5;
B = 4;

rho = lambda / (mu * m);

p_0 = p0(m, rho, B);
fprintf('p_0 : %.16e\n', p_0);

arr = zeros(1,B+1);
arr(1) = p_0;
for k = 1:1:B
    p_k = pk(m, rho, B, k, p_0);
    arr(k+1) = p_k;
    assignin('base', strcat('p_', num2str(k)), p_k);
    fprintf('%s : %.16e\n', strcat('p_', num2str(k)), p_k);
    % display(p_k);
end

somma = sum(arr);
fprintf('somma : %.16e\n', somma);

clear p_k;

